#*************        丁香园-NHANES     *************
#*************        Analyze Code      *************
#*************       Table 专项讲解     *************
#*************                          *************
#### 0.读取数据 & 准备环境 ####
library(gtsummary)
library(survey)
library(tidyverse)
setwd("F:/Rproject课程_5.3_Table1_专项讲解以及演示")

paper.data <- readRDS('Table1_Sample_Data.RData')
NHANES_design <- svydesign(data = paper.data, ids = ~SDMVPSU, strata = ~SDMVSTRA, nest = TRUE, weights = ~WTDRD1, survey.lonely.psu = "adjust") 

#### 1.Table1 基础版本，获得表格主体 ####

# by: 将 Table 进行分列、分队列呈现描述性统计的结果
# include: Table 中要呈现的变量，以及变量的排序
# label: 手动更改 Table 中变量的呈现内容，例如将 Age 这个变量在 Table 中呈现为 Age (years)
# type: 确定变量的类别
# statistic: 确定描述性变量呈现的内容，例如：均值、中位数等，
# 连续性变量常用的为：{median} ({p25}, {p75})，或者常用的："{mean} ({sd})"
# 离散性变量常用的为：{n_unweighted} ({p}%)

# digits: 修改 Table 中呈现的数值的小数位，可单独指定
# sort：改变变量取值排列的顺序，可以按照 "frequency" 或者 "alphanumeric" 进行排序, 或者是 everything() ~ "frequency"
# missing: "no", "ifany", "always"，是否在Table中呈现缺失值，默认取 no


tbl <- tbl_svysummary(NHANES_design,  
                      by = Sex,
                      include = c(Age, PIR, Race, BMI.group),
                      label = list(Age ~ 'Age (years)', PIR ~ 'PIR'),
                      # type = list(Age ~ "categorical"),
                      statistic = list(all_continuous()  ~ "{median} ({p25}, {p75})", 
                                       all_categorical() ~ "{n_unweighted} ({p}%)"),
                      digits = list(Age ~ 1, PIR ~ 2),  
                      sort = list(Race ~ "frequency", BMI.group ~ "alphanumeric" ),
                      missing = 'no') 
tbl



#### 2.Table1 进阶版本（一）：加上列：总队列的、P value 等信息 ####
# add_p(): 添加 p-values 列
# add_n(): 添加具有非缺失观察值总数的列，是否显示脚注，将 add_n 的关键参数进行呈现
# add_overall()：添加具有总体汇总统计信息的列
# add_q()：添加校正后的 p-values 结果，e.g. add_q(method = "bonferroni")
tbl_svysummary(NHANES_design,  
                      by = Sex,
                      include = c(Age, PIR, Race, BMI.group),
                      label = list(Age ~ 'Age (years)', PIR ~ 'PIR'),
                      # type = list(Age ~ "categorical"),
                      statistic = list(all_continuous()  ~ "{median} ({p25}, {p75})", 
                                       all_categorical() ~ "{n_unweighted} ({p}%)"),
                      digits = list(Age ~ 1, PIR ~ 2),  
                      sort = list(Race ~ "frequency", BMI.group ~ "alphanumeric" ),
                      missing = 'no') %>% 
  add_overall() %>%
  add_p()  %>%
  add_n(statistic = "{n}", # 默认为 "{n}, {N_miss_unweighted}
        # col_label = "**N**"
       # footnote = TRUE
       )
tbl


#### 3.Table1 进阶版本（二）：修改表的 Header 和脚注等 ####
# modify_spanning_header，stat_0 ~ NA，将 Overall 空出，至于这一列要叫什么名字，用show_header_names(tbl)获取
# modify_footnote：

tbl <- tbl_svysummary(NHANES_design,  
                      by = Sex,
                      include = c(Age, PIR, Race, BMI.group),
                      label = list(Age ~ 'Age (years)', PIR ~ 'PIR'),
                      # type = list(Age ~ "categorical"),
                      statistic = list(all_continuous()  ~ "{median} ({p25}, {p75})", 
                                       all_categorical() ~ "{n_unweighted} ({p}%)"),
                      digits = list(Age ~ 1, PIR ~ 2),  
                      sort = list(Race ~ "frequency", BMI.group ~ "alphanumeric" ),
                      missing = 'no') %>% 
  add_overall() %>%
  add_p() %>%
  add_n(statistic = "{N_nonmiss_unweighted}", # 默认为 "{n}, {N_miss_unweighted}
        col_label = "**N**",
        footnote = TRUE) %>% 
  modify_header(
    all_stat_cols() ~ "**{level}**, N = {n_unweighted} ({style_percent(p)}%)",
    p.value = "**P Value**") %>% # update the column header
  modify_spanning_header(
    stat_0 ~ NA, 
    update = all_stat_cols() ~ "**Sex**") %>% 
  modify_footnote(
    update = all_stat_cols() ~ 
      "median (IQR) for continuous; n (%) for categorical"
  )

show_header_names(tbl)
tbl

#### 4.Table1 进阶版本（三）：加粗/斜体设置 ####
# bold_labels(): 加粗变量
# italicize_levels(): 加斜体（嗯...有点奇怪，一般不要用，奇奇怪怪）
# bold_p(): 加粗 P value，可以设置 threshold，小于某个值的 p value 会变粗
tbl <- tbl_svysummary(NHANES_design,  
                      by = Sex,
                      include = c(Age, PIR, Race, BMI.group),
                      label = list(Age ~ 'Age (years)', PIR ~ 'PIR'),
                      # type = list(Age ~ "categorical"),
                      statistic = list(all_continuous()  ~ "{median} ({p25}, {p75})", 
                                       all_categorical() ~ "{n_unweighted} ({p}%)"),
                      digits = list(Age ~ 1, PIR ~ 2),  
                      sort = list(Race ~ "frequency", BMI.group ~ "alphanumeric" ),
                      missing = 'no') %>% 
  add_overall() %>%
  add_p() %>%
  add_n(statistic = "{N_nonmiss_unweighted}", # 默认为 "{n}, {N_miss_unweighted}
        col_label = "**N**",
        footnote = TRUE) %>% 
  modify_header(
    all_stat_cols() ~ "**{level}**, N = {n_unweighted} ({style_percent(p)}%)",
    p.value = "**P Value**") %>% # update the column header
  modify_spanning_header(
    stat_0 ~ NA, 
    update = all_stat_cols() ~ "**Sex**") %>% 
  modify_footnote(
    update = all_stat_cols() ~ "median (IQR) for continuous; n (%) for categorical") %>%
  bold_labels() %>%
  italicize_levels() %>%
  bold_p(0.05) 
tbl

#### 5.Table1 进阶版本（四）：左右叠加不同的 Table ####

tbl.all <- tbl_svysummary(NHANES_design,  
                          include = c(Age, PIR, Race, BMI.group),
                          label = list(Age ~ 'Age (years)', PIR ~ 'PIR'),
                          # type = list(Age ~ "categorical"),
                          statistic = list(all_continuous()  ~ "{median} ({p25}, {p75})", 
                                           all_categorical() ~ "{n_unweighted} ({p}%)"),
                          digits = list(Age ~ 1, PIR ~ 2),  
                          sort = list(Race ~ "frequency", BMI.group ~ "alphanumeric" ),
                          missing = 'no') %>%
  add_n(statistic = "{N_nonmiss_unweighted}", # 默认为 "{n}, {N_miss_unweighted}
        col_label = "**N**",
        footnote = TRUE) 


tbl.1 <- tbl_svysummary(NHANES_design,  
                        by = Sex,
                        include = c(Age, PIR, Race, BMI.group),
                        label = list(Age ~ 'Age (years)', PIR ~ 'PIR'),
                        # type = list(Age ~ "categorical"),
                        statistic = list(all_continuous()  ~ "{median} ({p25}, {p75})", 
                                         all_categorical() ~ "{n_unweighted} ({p}%)"),
                        digits = list(Age ~ 1, PIR ~ 2),  
                        sort = list(Race ~ "frequency", BMI.group ~ "alphanumeric" ),
                        missing = 'no') %>%
  add_p()


tbl.2 <- tbl_svysummary(NHANES_design,  
                        by = Age.group,
                        include = c(Age, PIR, Race, BMI.group),
                        label = list(Age ~ 'Age (years)', PIR ~ 'PIR'),
                        # type = list(Age ~ "categorical"),
                        statistic = list(all_continuous()  ~ "{median} ({p25}, {p75})", 
                                         all_categorical() ~ "{n_unweighted} ({p}%)"),
                        digits = list(Age ~ 1, PIR ~ 2),  
                        sort = list(Race ~ "frequency", BMI.group ~ "alphanumeric" ),
                        missing = 'no') %>%
  add_p() 

tbl_merge(
  tbls = list(tbl.all, tbl.1, tbl.2),
  tab_spanner = c("**Overall**", "**Sex Group**", "**Age Group**"))%>%
  modify_header(
    all_stat_cols() ~ "**{level}**, N = {n_unweighted} ({style_percent(p)}%)") %>% 
  bold_labels() 


